!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine resistivity(bx_resis, by_resis, bz_resis, en_resis, az_resis, ro, mx, my, mz, bx, by, bz, en,                           &
   x, y, dx, dy, nx, ny, max_et)
!==================================================================================================================================|

   implicit none

   integer(4), intent(in) :: nx, ny
   real(8), intent(in) :: x(nx), y(ny), dx, dy
   real(8), intent(inout) :: max_et
   real(8), dimension(nx, ny), intent(inout) :: bx_resis, by_resis, bz_resis, en_resis, az_resis, ro, mx, my, mz, bx, by, bz, en

   integer(4) :: i, j, ip1, im1, jp1, jm1

   real(8) :: dx2, dy2, dxdx, dydy
   real(8), dimension(nx, ny) :: jx, jy, jz, fxen_resis, fyen_resis, et

!----------------------------------------------------------------------------------------------------------------------------------|
!  initialization
!----------------------------------------------------------------------------------------------------------------------------------|
   dxdx = dx * dx
   dydy = dy * dy
   dx2 = dx * 2.0d0
   dy2 = dy * 2.0d0
   do j = 1, ny
   do i = 1, nx
      et(i, j) = 0.0d0
      jx(i, j) = 0.0d0
      jy(i, j) = 0.0d0
      jz(i, j) = 0.0d0
      bx_resis(i, j) = 0.0d0
      by_resis(i, j) = 0.0d0
      bz_resis(i, j) = 0.0d0
      en_resis(i, j) = 0.0d0
      az_resis(i, j) = 0.0d0
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  current density
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 2, ny - 1
      jp1 = j + 1
      jm1 = j - 1
      do i = 2, nx - 1
         ip1 = i + 1
         im1 = i - 1
         jx(i, j) = (bz(i, jp1) - bz(i, jm1)) / dy2
         jy(i, j) = - (bz(ip1, j) - bz(im1, j)) / dx2
         jz(i, j) = (by(ip1, j) - by(im1, j)) / dx2 - (bx(i, jp1) - bx(i, jm1)) / dy2
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  resistivity
!----------------------------------------------------------------------------------------------------------------------------------|
   call model_resis(et, ro, mx, my, mz, bx, by, bz, en, jx, jy, jz, x, y, nx, ny)

   max_et = max(maxval (et), max_et)

!----------------------------------------------------------------------------------------------------------------------------------|
!  variables update
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, ny
   do i = 1, nx
      fxen_resis(i, j) = et(i, j) * (by(i, j) * jz(i, j) - bz(i, j) * jy(i, j))
      fyen_resis(i, j) = et(i, j) * (bz(i, j) * jx(i, j) - bx(i, j) * jz(i, j))
   enddo
   enddo

   do j = 2, ny - 1
      jp1 = j + 1
      jm1 = j - 1
      do i = 2, nx - 1
         ip1 = i + 1
         im1 = i - 1
         bx_resis(i, j) = - (et(i, jp1) * jz(i, jp1) - et(i, jm1) * jz(i, jm1)) / dy2
         by_resis(i, j) = (et(ip1, j) * jz(ip1, j) - et(im1, j) * jz(im1, j)) / dx2
         bz_resis(i, j) = (et(i, jp1) * jx(i, jp1) - et(i, jm1) * jx(i, jm1)) / dy2 -                                              &
                          (et(ip1, j) * jy(ip1, j) - et(im1, j) * jy(im1, j)) / dx2
         en_resis(i, j) = (fxen_resis(ip1, j) - fxen_resis(im1, j)) / dx2 + (fyen_resis(i, jp1) - fyen_resis(i, jm1)) / dy2
         az_resis(i, j) = - et(i, j) * jz(i, j)
      enddo
   enddo

!   do j = 1, ny
!   do i = 1, nx
!      bx(i, j) = bx(i, j) + dt * bx_resis(i, j)
!      by(i, j) = by(i, j) + dt * by_resis(i, j)
!      bz(i, j) = bz(i, j) + dt * bz_resis(i, j)
!      en(i, j) = en(i, j) + dt * en_resis(i, j)
!      az(i, j) = az(i, j) + dt * az_resis(i, j)
!   enddo
!   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine resistivity
